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Abstract. The three-dimensional Thomas-Fermi (TF) model is used to simulate the variation of the d+d — > 
t + p cross-section at low impact energies, when the target deuterium nucleus is embedded in metallic or 
insulator environments. Comparison of the computational results to recent experiments demonstrates that 
even though the TF model can explain some increase in the low energy cross section for metallic host, a full 
explanation of the experimental results is still lacking. Possible reasons for the disagreement are discussed. 
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1 Introduction 

In a series of recent experiments [1-6] significant differ- 
ences had been found between the low-energy cross sec- 
tions of the d + d — > t + p reaction when the target nuclei 
are embedded in a metallic or an insulator environment. 
Apparently, the source of this effect is the conduction elec- 
tron distribution in the metallic lattice and its screening 
effects on the Coulomb barrier around the incident and 
target nuclei [7, 8] . Hints for similar behavior were found 
also in the electron capture rate of Be 7 [9, 10]. While in 
some of these experiments temperature dependence was 
also claimed [1,2], no such dependence was found in oth- 
ers [9,11]. 

This phenomenon can be briefly introduced as follows: 
The transparency, T(Ek), of a potential barrier to an in- 
cident particle having energy Ek is given by [12], 



T{E k ) = exp <^ -2 



2M r 



E k ] dr (1) 



where M is the reduced mass of the two particles, r/v - the 
nuclear radius, b - the classical turning point (CTP) and 
E p (r) — zeV(r) is the potential energy of the incident 
particle in the potential, V(r), generated by the target 
nucleus. Ze and ze are charges of the target and incident 
nuclei, respectively. 

If a pure Coulomb potential is substituted for the in- 
teraction potential, V(r) = Vc(r) = Ze / r, Eq. (1) repro- 
duces the Sommerfeld factor, 
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In (3) Vk = y / 2Ek/M is the velocity of the incident par- 
ticle. If the target nucleus is embedded in a solid target, 
the potential generated by nearby electrons, in addition 
to the bare Coulomb potential, has to be included as well, 



V(r) = V n (r)+V e (r) 



Ze 



K(r) 



(4) 



For low energy reactions the electronic part of the po- 
tential, V e (r), has a small but important contribution to 
the total potential. This contribution is amplified by the 
exponential factor of the transparency. 

In general, the electronic potential generated in a metal- 
lic lattice by the conduction electrons is a slowly varying 
function. Even for low energy particles the CTP of the 
projectile is so close to the target nucleus that one can 
fairly assume that along the integration path in (1) the 
electronic potential practically equals its value at the tar- 
get nucleus, V e (r) = 14(0). When this approximation is 
inserted back into (1) it gets the form, 



T{E k ) £* exp | -2 [ 
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This is equivalent to a reaction in a pure Coulomb 
potential with the kinetic energy, E^, replaced by Ek + U e , 
where 

U e = ze\V e (Q)\ (V e (0)<0) (6) 

The authors of Ref. [1] tried to explain the experi- 
mental results by means of the Debye statistical model, 
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which neglects the electrons degeneracy. A difficulty of 
their treatment is also the point that under the experi- 
mental conditions the Debye radius is significantly smaller 
than the atomic radius, therefore the basic assumptions of 
the model are not satisfied. 

The aim of the present paper is to compute the elec- 
trons spatial distribution in insulator and metallic envi- 
ronments for the d + d — > t+p reaction [Z = z = 1), using 
a three-dimensional Thomas-Fermi (TF) model, which ac- 
counts for the electrons' degeneracy, therefore presumably 
better describing the experimental conditions in a solid 
host of local high density conditions ("strongly-coupled 
plasma" ) . Using the TF model we try to infer the change 
in the transparency between metallic and insulating envi- 
ronments. To illustrate our method, we shall focus on the 
case of deuterium embedded in a copper lattice for which 
experimental results are available [1]. 



2 The model 

2.1 Basic data 

In the experiment, Ref. [1], the deuterium atoms consist 
only a small part of the target - about 11 copper atoms 
for each deuterium atom [1] - one can, therefore, safely 
assume that the presence of the deuterium atoms does 
not significantly modify the copper lattice properties. 

Our first step is to find the volume available for the 
deuterium atoms. This is carried out by means of the 
QEOS method [13], which is frequently used to find the 
equation-of-state of various materials, and was found to 
give accurate results. For our purposes its main advan- 
tage is that QEOS can provide the volume per atom sep- 
arately for each component in a mixture of materials. As- 
suming a deuterium/copper solid material with 9% deu- 
terium and 91% copper at 8.93 g/cm 3 (solid copper spe- 
cific gravity), QEOS predicts that the volumes of the deu- 
terium (copper) atoms in the target are VoId (Volcu) = 
1.79 (11.7)10 -24 cm 3 / atom, i.e. , the average volume avail- 
able for a copper atom is by a factor of ~ 6.5 larger than 
that of a deuterium atom. If these atomic volumes are as- 
sumed to have the shape of a spherical enclosure, called in 
the following the Ion Sphere (IS), then the corresponding 
IS radii are R l:D {R lfi u) = 0.717(1.36) 1CT 8 cm. 

The atomic structure of copper is [Ar] 3d 10 4s. In a 
pure copper metal lattice there are at most 11 electrons 
per atom in the conduction band. Their average density 
is n e = 11/11.710" 24 = 9.4T0 23 electrons / cm 3 . If the 
small amount of deuterium atoms does not significantly 
change this distribution, then the same density prevails 
within the deuterium IS as well, generating n e Volo = 
9.4T0 23 x 1.79T0 -24 = 1.7 extra electrons inside the deu- 
terium IS. Thus, together with its own electron, the deu- 
terium IS contains, on the average, N e = 2.7 electrons. 
Changing the Cu/D target density by ±10% changes this 
quantity only by ±5%. In our computations we have used 
N e , the number of electrons in the deuterium IS, as an 
adjustable parameter. 



2.2 The modelling of the electronic potential in 
insulating and metallic environments 

The electronic potential of a free deuterium atom in high 
density insulating environment was calculated by the com- 
puter program RELDIR, which solves the relativistic Dirac 
equation of a deuterium (or any other) atom in a finite ra- 
dius Ion Sphere with any number of free electrons within 
the IS. The results of this computation are the bound and 
free electrons wavefunctions, their spatial density, and the 
corresponding potentials. 

As a first step, the computation was carried out for 
a solid density deuterium lattice, which is known to be 
an insulator, (0.202g/cm 3 , 6. 08T0 22 atoms / cm 3 , Ri = 
1.58T0 _8 cm). From the result of RELDIR for the elec- 
tronic potential one obtains U e = 20 eV, in accordance 
with the experimental results of Ref. [1] in an insulating 
environment. This result was used as the basic reference 
for comparison. 

The evaluation of the electronic potential near a deu- 
terium atom embedded in a metallic lattice was carried 
out by means of the Thomas-Fermi (TF) model in con- 
junction with the Born-Oppenheimer (BO) approxima- 
tion. Owing to the high density environment, a Fermi- 
Dirac statistics for the electrons, as used by the TF model, 
provides a better approximation for the electrons spatial 
and energy distributions inside the copper lattice than the 
Debye model [1,2] which is more appropriate for weakly- 
coupled low-density matter. In this context, the physical 
meaning of the BO approximation is that nearby electrons 
rapidly adjust their local distribution to any change in the 
nuclei positions, so that the electron distribution depends 
only on the instantaneous distance between the two nuclei, 
but not on their relative motion. 

The TF model for the electrons consists of a set of 
three equations. The first one is the Poisson equation, 



vV e (r) = 47ren e (e F ;r) 



(7) 



whose solution is the electronic potential, V e (r), when the 
electrons spatial density, n e (e F ',r), is known. In (7) e F is 
the Fermi energy. The second equation is the Fermi-Dirac 
distribution of the electrons inside the IS, 



n e (e F ;r) 



xfi/2 



1 

2^ 



s F + eV{r) _ 
kT 
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3/2 



3V(V) 
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(8) 



which provides the electron density as function of the total 
potential V(r) = V^r) + V e (r). In (8), kT is the temper- 
ature of the target material (= 300 K) in energy units, 
and, 

f°° v 1/2 dv 

is the incomplete Fermi-Dirac integral [15]. The boundary 
condition applied to Eqs. (7) and (8) is Gauss' theorem, 
//s surface E ' dS = 4ire N e . Finally, the Fermi energy, e F , 
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is computed from the condition that the total charge inside 
the IS equals N e , 

N e = [ n e (e F ;r)d 3 r (10) 

J is volume 

Simultaneous solution of Eqs. (8)-(10) yields the elec- 
tronic potential, V e (r), the electron density, n e (sF\ r), and 
the Fermi energy, sf- When these solutions are inserted 
back into (1) and (4), one finds the CTP, b(N e ,E k ), and 
the transparency, T(N e , as function of the number of 
electrons inside the IS and the incident particle's kinetic 
energy. 



The internal iterations calculated the electron density and po- 
tential through equations (12)-(14) and (8)-(10). When these 
have been found, a second iterative procedure was applied to 
solve b from the condition, 

E p (b) = E k (15) 

This double iterative method was continued until conver- 
gence was achieved for all the parameters. 



3 Results 

3.1 The electron density 



2.3 Computational details 

The solution of (7) can be rewritten as [14], 



V e (r,6,if) = -e 



is volume 

Ri 

r 12 d 



(11) 



6' d6' 



r v w 

J^o k — r'j 



Obviously, the electron distribution has cylindrical sym- 
metry around the line connecting the two nuclei, therefore 
n e (r',9') is independent of ip' . Moreover, the charge dis- 
tribution has also a reflection symmetry around the plane 
perpendicular to this line at halfway between the nuclei. 
Using these symmetry properties, the integration over ip' 
can be carried out analytically, thereby reducing the triple 
integral in (11) to a double one, see details in Ref. [14]. 
Finally, transforming the coordinate system center onto 
the target nucleus, Eq. (11) gets the form, 



V e (R,fi) = -4e / R' 2 dR' 



>, f R ' r' 2 dR' r 



dp! n e (R' , fi' 
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Figure 1 shows the density of the electrons around the incident 
and the target nuclei, when the center-of-mass energy of the 
two nuclei is 4000 eV. The figure clearly shows the polarization 
of the electrons near the nuclei. We recall that the TF model 
predicts that the electron density diverges as r~ 3 ^ 2 near the 
nucleus [15]. 
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Fig. 1. The distribution of the electrons around the two nuclei. 



A± = [(b/2 + Rcos8)± (b/2 + R' cos6')] 2 

+ (RsmO) 2 + (R' sine')' 2 (13) 

B = 2R R' sin 6 sin 9' (14) 

In (12) R is the radial distance of the field point from the 
target nucleus, fi = cos#, K(x) is the Jacobi elliptic integral 
[16], and /i m i n = max(—b/2R', —1). The reduction of the 
triple integral in (11) to the shape of (12) reduced the compu- 
tational resources and improved the accuracy of the numerical 
procedure. 

The numerical process was greatly complicated by the fact 
that the CTP, b, not only depends on the electron density, 
n e (eF',r), but also determines it. It was, therefore, necessary 
to carry out a doubly iterative method for the computations. 



The TF model predicts an accumulation of the electrons 
along the line connecting the nuclei, with a saddle point at 
halfway between them, a fact which enhances the screening 
effect. 



3.2 The electronic potential 

Figure 2 shows the spatial distribution of the electronic poten- 
tial in the same region as in Fig. 1. In the interesting part of 
the field, i.e. in the space around the two nuclei, the poten- 
tial has very slow variation, which justifies the approximation 
V e (r) = F e (0), to an accuracy of a few percents. 

The value of the potential on the nuclei, U e = |eI4(0)|, is 
displayed on Fig. 3. This is, in fact, the quantity that is mea- 
sured experimentally. The figure shows that U e fits excellently 
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Fig. 2. The distribution of the electronic potential, V e (r), 
around the two nuclei. 
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Fig. 3. The value of the electronic potential, V e (r), on the 
target nucleus. 

a linear behavior as function of N e , with a small contribution 
from the incident energy, 

3D : U e = |eV e (0)| = 27.2 iV e + 80.3 + 3.ri0" 4 £ fe eV (16) 

The last term is always less than 4% in the range of interest. 
In order to compare the importance of a 3D modelling, we 
have developed also a one-dimensional TF program that solves 
the same problem in a radially symmetrical environment. The 
results of this program is also plotted on Fig. 3. The TF ID 
code yielded, 

ID : U e — |eV e (0)j = 27.2 7V e + 33.6 eV (17) 

The difference between these two results originates from 
the fact that the electrons accumulation between the nuclei in 
a 3D model is larger than in the ID case. This difference clearly 



indicates the importance of a three-dimensional treatment of 
the problem. 

For N e = 3, E k = 4000 eV Eq. (16) gives U e = 163 eV. 
This has to be compared to the experimental result, U e = 
470 ± 50 eV , which is by a factor of 2.9 ± 0.3 larger than the 
result of the TF 3D model. It has also to be compared to 
U e = 20 eV as computed from RELDIR for the insulator case, 
see above. 

3.3 The Fermi energy 

The Fermi energy, ef, strongly depends on the number of elec- 
trons in the IS, N e , and is, of course, independent on the incom- 
ing particle's energy. The behavior of ef vs. N e is illustrated 
in Fig. 4. Within the range of our computations ef turned out 
to be a linear function of N e , 

3D: ef = -43.7 + 38.8 N e (eV) (18) 

with a hint for a slight quadratic curvature. For the ID case 
we have found a significantly higher result, 

ID: ef = -15.9 + 35.3 N e (eV) (19) 

It should be emphasized, that the Fermi energy is positive 
(except for the case of N B = 1, namely, the case of an isolated 
deuterium atom), and is much higher than the target temper- 
ature (e f > kT = 300 if = 0.025 eV), for all the cases. This 
means that at the temperatures of the experiments (room tem- 
perature and below) there is no reason for measurable temper- 
ature variations. In fact, we have repeated our computations 
with T = 3000 K , but this order of magnitude change in the 
temperature modified the results by less than 0.5% - as ex- 
pected. 

3.4 The classical turning point 

The CTP fits, with excellent accuracy, a function of the form, 
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Fig. 4. Fermi energy as function of the number of electrons in- 
side the ion sphere, for 3D and ID Thomas-Fermi simulations. 
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where e 2 / do = 27.211 eV, (ao is the Bohr radius). This has 
exactly the form of the Coulomb CTP, with the projectile's 
kinetic energy modified according to Eq. (5). Fitting the com- 
putational results to the form of (20) provides, 



U e = 27.7 N e + 76.4 + 5.76' W~ 4 E k eV 



(21) 



The agreement of this U e with the values obtained from the 
electronic potential, Eqs. (16), can be regarded as very good. 
The difference between the two results stems from the fact 
that |eV e (r)| is not exactly constant along the line connecting 
the two nuclei. And again, this result is in discrepancy with 
the experimental value, but is much higher than the insulator 



3.5 The transparency 

At low impact energies the screened astrophysical factor, S(E), 
is enhanced relative to the unscreened one by T(E)/Tc(E) [1], 
see Eqs. (1) and (2). Fig. 5 shows our results for this factor 
as function of the center of mass energy for the d + d — > t + p 
reaction, when the target deuterium is embedded in a copper 
lattice. The experimental results of Ref. [1] are also displayed 
in Fig. 5. The dashed curve represents the bare S(E) factor [1]. 
While we note that the value N e = 3 is a reasonable estimate 
for the conducting electrons in the deuterium IS embedded in 
Cu, even for N e — 8 the 3-dimcnsional TF model underesti- 
mates the experimental low-energy increase of the astrophysi- 
cal factor. The low U e predicted by the model relative to the 
experiment is another manifestation of the same fact. In fact, 
one needs 11 conduction electrons inside the deuterium IS to 
get agreement with the experimental results. 

In the range of the energies used in our computations U e < 
0.1-E. Denoting £ = U e / E, £ can be assumed to be a small 
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Fig. 5. Comparison between the results of the 3D TF model 
for the astrophysical factor, S(E), of the reaction d + d — ► t+p, 
when the target nuclei are embedded in a copper substrate, and 
the experimental results. The experimental points are taken 
from [1]. The dotted line presented the theoretically extrapo- 
lated "bare nucleon" S(E) factor 



quantity. Using first order expansion, an analytical estimate 
can be developed for the change in the value of T(E) / Tc{E) 
caused by a change b — » 6/ (1 + £) in the CTP, see Eq. (20). 
The ratio of the screened to Coulomb transparencies becomes, 



T(E)/T C = exp{5G} 



where, 



(22) 



(23) 



For E = 4000 eV, N e = 3 and U e = 163 eV this formula 
predicts, 

T{E)/T c = exp{6G} = lA5 (24) 
in contrast to the experimental result of T(E) / Tc « 2.0±0.2 . 



4 Discussion 

In this paper we present computational results from a three- 
dimensional TF model about the modification of low-energy 
cross-section for the d + d — > t + p reaction, when the target 
nucleus is embedded in a copper lattice. Our computational 
result of U e = 163 eV is lower by a factor of ~3 from the ex- 
perimental result of U e = 470 ± 50 eV, but still is substantially 
higher than the cross-section in an insulator (solid deuterium 
lattice), Ue = 20 eV. 

In order to see whether such a difference holds true for 
other target lattices as well, in Fig. 6 we plotted the results 
of all the experiments in metal lattices published in Ref. [1] 
in comparison to ours. The 3D TF results in the figure may 
be shifted ±5% up or down, due to differences in the local IS 
volume of the deuterium in the various lattices. Figure 6 ex- 
hibits a consistent difference between the computational and 
the experimental results. As the Thomas-Fermi is a highly suc- 
cessful model, well fitted to high density matter, this insistent 
difference is of some surprise. 
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Fig. 6. The experimental values of U e , drawn from the results 
of [1] , compared to the predictions of the 3D TF model. 
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Obviously, the central reason for the disagreement is the 
electronic potential on the target nucleus, see Eq. (6). In our 
opinion, the semiclassical nature of the TF model cannot be 
the reason for this difference, because the quantum mechanical 
behavior of the conduction electrons seems to be unimportant 
in the present problem. 

On the other hand, the TF model assumption of a perfect 
spherical symmetry, and the target nucleus location at the cen- 
ter of a well-defined ion sphere, may oversimplify the real sit- 
uation. In the experiment, the deuterium nuclei are implanted 
into the metal lattice by bombardment of 10 keV deuterons 
into the metal foil [2] . It is well known that this technique does 
not necessarily deposit the stopped deuterium into a spheri- 
cally symmetrical environment. 

Another contribution to the discrepancy may be the local 
conduction electron density fluctuations around the target nu- 
cleus. As the transparency depends exponentially on the local 
electron potential, a small number of nuclei, with an instan- 
taneous large electron screening, has greater influence on the 
cross-section than the other nuclei with the average electronic 
screening. 

Finally, it is possible that QEOS underestimates the deu- 
terium IS volume under the specific experimental conditions. 
Larger IS volume, with more conduction electrons, may have 
better agreement with the experiment. 
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